Produced with R version 4.0.0.

Preparations

Apply TSIR model to the individual towns.

Susceptible reconstruction.

Smoothing spline regression of cumulative cases on cumulative births to estimate under-reporting and residuals. Correct for under-reporting: I = cases/ur.

Regress cumulative cases on cumulative births to estimate under-reporting (ur) and susceptible depletion (z) I is estimated as cases/ur.

Now reconstruct the susceptibles (via the residuals z).

Does the preceding scale the z variable appropriately? I.e., does it take under-reporting into account?

Fit TSIR

First, set up the data matrix for the regression. This involves lagging the I and z variables. Population size is assumed constant at its median value.

We’ll start by estimating the mean susceptible fraction, assuming a single global value for this parameter. Some towns have very high changes in susceptible fraction; exclude these.

## [1] 233

Now profile on the fraction, \(\sigma\), of susceptibles in the population. This assumes a single, global value of \(\sigma\). NOTE: the \(\beta\) are here scaled by population size.

Our global sigma estimate is \(\sigma=0.0512\). Now, we assume this value of \(\sigma\) and fit TSIR to each of the towns individually using linear regression.

Just for interest, let’s plot \(R_0\) as a function of city size.

  • Shat is fitted to each individually.
  • Sbar is fitted globally (with some towns excluded).

Coupling the cities.

Let \(X_{i}(t)\) be the observation at time \(t\) in city \(i\). We assume that \[X_{i}(t+1) \sim \mathrm{Poisson}\left(\lambda_i(t)\right)\] or, alternatively, \[X_{i}(t+1) \sim \mathrm{NegBinom}\left(\lambda_i(t),\frac{1}{\psi}\right),\] where \(\psi\) is an overdispersion parameter. In the latter, we have parameterized the negative binomial distribution so that \(\mathbb{E}\left[{X_i(t+1)}\right]=\lambda_i(t)\) and \(\mathrm{Var}\left[{X_i(t+1)}\right]=\lambda_i(t)+\psi\,\lambda_i(t)^2\).

When \(I_i(t)=0\), we have that \[\lambda_i(t) = \beta_i(t)\,S_i(t)\,\iota_i(t)^\alpha.\] In the above, \(\beta\) is constructed by fitting the TSIR model to each city independently and the susceptible pool, \(S\), is reconstructed using TSIR methods. The quantity \(\iota\) is the import rate, which we estimate using a variety of different models.

The following computes \(y\), \(S\), \(\beta\), and the matrix of reciprocal distances. It also picks out the relevant observations. These are the ones for which the preceding week saw zero cases.

Gravity model

The gravity model (with intercept) is \[\iota_i=N_i^{\tau_1} \left(\theta\,\sum_{j\ne i}\! N_j^{\tau_2} d_{ij}^{-\rho}\,\frac{I_j}{N_j}+\phi\right).\] Let \[Q_{ij}=\begin{cases}N_i^{\tau_2}\,d_{ij}^{-\rho}, &i\ne j\\0, &i=j\end{cases}\] and \[y_{i}=\frac{I_i}{N_i}.\]

Expressed compactly, the gravity model is \[\iota = \theta\,\mathrm{diag}(N^{\tau_1})\,Q^T\,y+\phi\,N^{\tau_1}.\]

We use R’s element-recycling feature, which allows us to write \[\mathrm{diag}(A)\,B=A\;*\;B.\]

The negative log likelihood function for the gravity model is:

Now we’ll do a likelihood profile over \(\tau_1\) and \(\tau_2\).

bake(file="gravity.rds",{
  
  grd <- profileDesign(
    tau1=seq(0.1, 1.4, length=25),
    tau2=seq(-1, 1, length=25),
    lower=c(psi=log(5),theta=log(1e-8),phi=log(1e-6),rho=log(0.9)),
    upper=c(psi=log(7),theta=log(1),phi=log(1e-2),rho=log(1.7)),
    nprof=10
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau1","tau2")
      formals(likfnGravity) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnGravity,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  results %>%
    extract(sapply(results,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  results %>%
    extract(!sapply(results,inherits,"try-error")) %>%
    ldply() %>%
    ddply(~tau1+tau2,subset,loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau1","tau2")
      formals(likfnGravity) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnGravity,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  
  attr(results,"nproc") <- getDoParWorkers()
  results
}) -> results

The above was completed in 49.1 min on a 250-core cluster. Now we refine to obtain the global MLE.

## list()

Regarding NLOPT return values:

Successful termination (positive return values)

  • NLOPT_SUCCESS = 1 Generic success return value.
  • NLOPT_STOPVAL_REACHED = 2 Optimization stopped because stopval (above) was reached.
  • NLOPT_FTOL_REACHED = 3 Optimization stopped because ftol_rel or ftol_abs (above) was reached.
  • NLOPT_XTOL_REACHED = 4 Optimization stopped because xtol_rel or xtol_abs (above) was reached.
  • NLOPT_MAXEVAL_REACHED = 5 Optimization stopped because maxeval (above) was reached.
  • NLOPT_MAXTIME_REACHED = 6 Optimization stopped because maxtime (above) was reached.

Error codes (negative return values)

  • NLOPT_FAILURE = -1 Generic failure code.
  • NLOPT_INVALID_ARGS = -2 Invalid arguments (e.g. lower bounds are bigger than upper bounds, an unknown algorithm was specified, etcetera).
  • NLOPT_OUT_OF_MEMORY = -3 Ran out of memory.
  • NLOPT_ROUNDOFF_LIMITED = -4 Halted because roundoff errors limited progress. (In this case, the optimization still typically returns a useful result.)
  • NLOPT_FORCED_STOP = -5 Halted because of a forced termination: the user called nlopt_force_stop(opt) on the optimization’s nlopt_opt object opt from the user’s objective function or constraints.

We look for evidence of overdispersion by computing the scale parameter for the negative binomial model.

Xia model

The model of Xia et al. (2004) is \[\iota_i=N_i^{\tau_1}\,\left(\theta\,\sum_{j\ne i}\! I_j^{\tau_2} d_{ij}^{-\rho} + \phi\right).\] Let \[Q_{ij}=\begin{cases}N_i^{\tau_2}\,d_{ij}^{-\rho}, &i\ne j\\0, &i=j\end{cases}\] and \[y_{i}=\frac{I_i}{N_i}.\]

Expressed compactly, the Xia et al. (2004) model is \[\iota = \theta\,\mathrm{diag}(N^{\tau_1})\,Q^T\,y^{\tau_2}+\phi\,N^{\tau_1}.\]

The negative log likelihood function for the Xia et al. (2004) model is:

Again, a profile computation.

The above was completed in 10.6 min on a 250-core cluster.

## list()

Mean field model

By setting \(\rho = 0\) and \(\tau_1=\tau_2=1\) in the gravity model, we obtain mean-field model.

The negative log likelihood function for the mean-field model is:

Now we find the MLE.

bake(file="meanfield.rds",{

  grd <- sobolDesign(
    lower=c(psi=log(5),theta=log(1e-8),phi=log(1e-6)),
    upper=c(psi=log(7),theta=log(1),phi=log(1e-2)),
    nseq=250
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c()
      formals(likfnMeanField) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnMeanField,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  results %>%
    extract(sapply(results,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  results %>%
    extract(!sapply(results,inherits,"try-error")) %>%
    ldply() %>%
    subset(loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c()
      formals(likfnMeanField) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnMeanField,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  attr(results,"nproc") <- getDoParWorkers()
  results
  
}) -> results

The above was completed in 2.4 min on a 250-core cluster. Now we refine to obtain the global MLE.

## list()

Diffusion model

By setting \(\tau_1 = \tau_2 = 0\) in the gravity model, we obtain a diffusion model, which couples cities in a way that depends only on the distance between them.

The negative log likelihood function for the diffusion model is:

Now we’ll do a likelihood profile over \(\tau_1\) and \(\tau_2\).

bake(file="diffusion.rds",{

  grd <- profileDesign(
    rho = seq(log(1),log(3),length=50),
    lower=c(psi=log(5),theta=log(1e-8),phi=log(1e-6)),
    upper=c(psi=log(7),theta=log(1),phi=log(1e-2)),
    nprof=10
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("rho")
      formals(likfnDiffusion) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnDiffusion,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  results %>%
    extract(sapply(results,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  results %>%
    extract(!sapply(results,inherits,"try-error")) %>%
    ldply() %>%
    ddply(~rho,subset,loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("rho")
      formals(likfnDiffusion) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnDiffusion,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  attr(results,"nproc") <- getDoParWorkers()
  results
  
}) -> results

The above was completed in 4.4 min on a 250-core cluster. Now we refine to obtain the global MLE.

## list()

Competing destinations model

The competing destinations model is \[\iota_i=N_i^{\tau_1}\,\left(\theta\,\sum_j\frac{N_j^{\tau_2}}{d_{ij}^\rho}\,\left(\sum_{k \ne i, j}\frac{N_k^{\tau_2}}{d_{jk}^\rho}\right)^\delta\,\frac{I_j}{N_j}+\phi\right).\]

Let \[Q_{ij}=\begin{cases}{N_i^{\tau_2}}{d_{ij}^{-\rho}}, &i\ne j\\0, &i=j\end{cases}\] and \[R_{ji}=\sum_{k \ne i, j}{N_k^{\tau_2}}{d_{jk}^{-\rho}}=\sum_{k\ne i,j}Q_{kj}=\sum_{k\ne i}Q_{kj}=\sum_{k}Q_{kj}-Q_{ij}.\] This implies \[R^T=(\mathbb{1}\,\mathbb{1}^T-I)\,Q \qquad \Longleftrightarrow \qquad R=Q^T\,(\mathbb{1}\,\mathbb{1}^T-I)\] and \[\iota_i=N_i^{\tau_1}\,\left(\theta\,\sum_j Q_{ji}\,R_{ji}^\delta\,y_{j}+\phi\right).\]

The negative log likelihood function for the competing destinations model is:

We compute the profile likelihood as before.

bake(file="compdest.rds",{
  
  grd <- profileDesign(
    tau1=seq(0.1, 1.4, length=25),
    tau2=seq(-1, 1, length=25),
    lower=c(psi=log(5),theta=log(0.1),phi=log(1e-6),rho=log(1.5),delta=-3),
    upper=c(psi=log(7),theta=log(45),phi=log(1e-2),rho=log(30),delta=0),
    nprof=20
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau1","tau2")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  results %>%
    extract(sapply(results,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  results %>%
    extract(!sapply(results,inherits,"try-error")) %>%
    ldply() %>%
    ddply(~tau1+tau2,subset,loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau1","tau2")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> results
  
  attr(results,"nproc") <- getDoParWorkers()
  
  results
}) -> results

The above was completed in 171.3 min on a 250-core cluster.

## list()

One-D profiles

\(\delta\)

bake(file="compdest_delta.rds",{
  
  grd <- profileDesign(
    delta=seq(-1.6,0,length=250),
    lower=c(tau1=0.1,tau2=-1,psi=log(5),theta=log(0.1),phi=log(1e-6),rho=log(1.5)),
    upper=c(tau1=1.4,tau2=1,psi=log(7),theta=log(45),phi=log(1e-2),rho=log(30)),
    nprof=50
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("delta")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> res
  
  res %>%
    extract(sapply(res,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  res %>%
    extract(!sapply(res,inherits,"try-error")) %>%
    ldply() %>%
    ddply(~delta,subset,loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("delta")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> res1
  
  attr(res1,"nproc") <- getDoParWorkers()
  
  res1
}) -> res1

The above was completed in 260 min on a 250-core cluster.

\(\tau_1\)

bake(file="compdest_tau1.rds",{
  
  grd <- profileDesign(
    tau1=seq(0.1,1.4,length=250),
    lower=c(tau2=-1,psi=log(5),theta=log(0.1),phi=log(1e-6),rho=log(1.5),delta=-3),
    upper=c(tau2=1,psi=log(7),theta=log(45),phi=log(1e-2),rho=log(30),delta=0),
    nprof=50
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau1")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> res
  
  res %>%
    extract(sapply(res,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  res %>%
    extract(!sapply(res,inherits,"try-error")) %>%
    ldply() %>%
    ddply(~tau1,subset,loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau1")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> res2
  
  attr(res2,"nproc") <- getDoParWorkers()
  
  res2
}) -> res2

The above was completed in 235.6 min on a 250-core cluster.

\(\tau_2\)

bake(file="compdest_tau2.rds",{
  
  grd <- profileDesign(
    tau2=seq(-1,1,length=250),
    lower=c(tau1=0.1,psi=log(5),theta=log(0.1),phi=log(1e-6),rho=log(1.5),delta=-3),
    upper=c(tau1=1.4,psi=log(7),theta=log(45),phi=log(1e-2),rho=log(30),delta=0),
    nprof=50
  )
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau2")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1000,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> res
  
  res %>%
    extract(sapply(res,inherits,"try-error")) %>%
    sapply(as.character) %>%
    unique() %>%
    print()
  
  res %>%
    extract(!sapply(res,inherits,"try-error")) %>%
    ldply() %>%
    ddply(~tau2,subset,loglik==max(loglik,na.rm=TRUE)) -> grd
  
  foreach(start=iter(grd,"row"),
    .errorhandling="pass",
    .inorder=FALSE,
    .packages=c("nloptr","magrittr")
  ) %dopar% try(
    {
      fixed <- c("tau2")
      formals(likfnCompDest) %>% names() %>% setdiff(fixed) -> est
      nloptr(unlist(start[est]),
        function(x){
          do.call(likfnCompDest,
            c(start[fixed],setNames(as.list(x),est)))
        },
        opts=list(
          algorithm="NLOPT_LN_SBPLX",
          ftol_abs=1e-5,
          xtol_rel=1e-7,
          maxeval=10000)
      ) -> fit
      fit$solution %>% setNames(est) %>%
        c(unlist(start[fixed]),loglik=-fit$objective) %>%
        as.list() %>% as.data.frame() %>%
        cbind(conv=fit$status) -> res
    }
  ) -> res3
  
  attr(res3,"nproc") <- getDoParWorkers()
  
  res3
}) -> res3

The above was completed in 279.8 min on a 250-core cluster.

Stouffer model

Let \(i\) be the recipient town; \(j\), the donor. Let \(S(i,j)\) be the collection of towns closer to town \(i\) than \(j\) is. That is, \(S(i,j) = \{k: k\ne i \ \&\ d(i,k) \le d(i,j)\}\).

\[\iota_i = N_i^{\tau_1}\,\left(\theta\,\sum_j\!\left(\frac{N_j}{\sum_{k\in S(i,j)}{N_k}}\right)^{\tau_2}\,\frac{I_j}{N_j}+\phi\right).\]

The negative log likelihood function for the Xia et al. (2004) model is:

We compute the profile likelihood as before.

The above was completed in 4.5 min on a 250-core cluster.

## list()

Variant: Stouffer model with recipient included

Let \(i\) be the recipient town; \(j\), the donor. Let \(S'(i,j)\) be the collection of towns closer to town \(i\) than \(j\) is. That is, \(S'(i,j) = \{k: d(i,k) \le d(i,j)\}\). Note that \(S'(i,j)\) includes \(i\), whilst \(S(i,j)\) does not.

\[\iota_i = N_i^{\tau_1}\,\left(\theta\,\sum_j\!\left(\frac{N_j}{\sum_{k\in S'(i,j)}{N_k}}\right)^{\tau_2}\,\frac{I_j}{N_j}+\phi\right)\]

The profile likelihood computation.

The above was completed in 4.5 min on a 250-core cluster.

## list()

Radiation model

\[\iota_i=\theta \sum_j N_j \frac{N_j N_i}{(N_j + \sum_{k \in S(i,j)} N_k)(N_j + N_i + \sum_{k \in S(i,j)} N_k)} \frac{I_j}{N_j}\]

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = likfnRadiation, start = list(theta = log(1), 
##     psi = log(5)), method = "Nelder-Mead", skip.hessian = FALSE, 
##     control = list(trace = 0, maxit = 10000))
## 
## Coefficients:
##         Estimate Std. Error z value     Pr(z)    
## theta -2.0846191  0.0069030 -301.99 < 2.2e-16 ***
## psi    1.7421955  0.0071838  242.52 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 400456.9
## 'log Lik.' -200228.5 (df=2)

Variant: radiation model with recipient included

\[\iota_i=\theta \sum_j N_j \frac{N_j N_i}{(N_j + \sum_{k \in S'(i,j)} N_k)(N_j + N_i + \sum_{k \in S'(i,j)} N_k)} \frac{I_j}{N_j}\]

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = likfnRadiation1, start = list(theta = log(1), 
##     psi = log(5)), method = "Nelder-Mead", skip.hessian = FALSE, 
##     control = list(trace = 0, maxit = 10000))
## 
## Coefficients:
##         Estimate Std. Error z value     Pr(z)    
## theta -1.8844494  0.0067635 -278.62 < 2.2e-16 ***
## psi    1.7252532  0.0072161  239.08 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 399364.3
## 'log Lik.' -199682.2 (df=2)

Extended radiation model

We extend the radiation variant model by allowing a background importation rate proportional to city size.

\[\iota_i=\theta \sum_j N_j \frac{N_j N_i}{(N_j + \sum_{k \in S'(i,j)} N_k)(N_j + N_i + \sum_{k \in S'(i,j)} N_k)} \frac{I_j}{N_j}+\phi\,N_i\]

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = likfnXRad1, start = list(theta = log(1), psi = log(5), 
##     phi = log(1e-05)), method = "Nelder-Mead", skip.hessian = FALSE, 
##     control = list(trace = 0, maxit = 10000))
## 
## Coefficients:
##          Estimate  Std. Error z value     Pr(z)    
## theta  -2.4559288   0.0111829 -219.62 < 2.2e-16 ***
## psi     1.6285242   0.0073544  221.43 < 2.2e-16 ***
## phi   -11.6821427   0.0183598 -636.29 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 393261.7
## 'log Lik.' -196630.9 (df=3)

Model comparison

model \(\ell\) \(\mathrm{QAIC}\) \(\hat{c}\) \(\Delta\!\mathrm{QAIC}\) \(\log_{10}{\theta}\) \(\log_{10}{\phi}\) \(\rho\) \(\tau_1\) \(\tau_2\) \(\delta\) \(\psi\)
SV -196466.8 137786.2 2.852 0.000 -0.796 -4.563 0.883 1.730 5.068
XR -196630.9 137897.2 3.186 111.021 -1.067 -5.073 5.096
CD -196654.3 137921.7 2.887 135.498 0.439 -3.699 1.860 0.661 0.501 -1.006 5.077
S -196778.9 138005.0 2.764 218.836 -0.820 -4.291 0.819 1.440 5.140
G -197733.0 138676.1 2.876 889.920 -3.453 -3.445 1.516 0.614 0.134 5.265
X -198028.0 138883.0 2.851 1096.833 -6.023 -6.303 0.820 0.607 0.433 5.327
RV -199682.2 140035.0 3.953 2248.813 -0.818 5.614
R -200228.5 140418.1 4.372 2631.918 -0.905 5.710
MF -200369.7 140523.2 4.098 2736.995 -8.769 -5.083 5.745
D -202735.5 142180.2 2.023 4394.008 -0.719 -0.788 1.942 6.263

Diagnostics

Import and export rates by population size

We can predict importation rates from fitted models

Likelihood differences

dat %>% acast(town~year+biweek,value.var="cases") -> obs1

likGravity <- function (mle) {
  npar <- 6
  mle %$% {
    theta <- exp(theta)
    phi <- exp(phi)
    rho <- exp(rho)
    Q <- (N^tau2)*(dd^rho)
    iota <- (N^tau1)*(theta*crossprod(Q,ylag)+phi)
    lambda <- betaS*iota[relevant]^alpha
    ll <- dnbinom(x=obs,mu=lambda,size=exp(-psi),log=TRUE)
    rv <- array(0,dim=dim(iota),dimnames=dimnames(iota))
    rv[relevant] <- ll
    rowSums(rv)/rowSums(rv!=0)
  }
}

likCompDest <- function (mle) {
  mle %$% {
    theta <- exp(theta)
    phi <- exp(phi)
    rho <- exp(rho)
    Q <- (N^tau2)*(dd^rho)
    R <- crossprod(Q,iii)^delta
    iota <- (N^tau1)*(theta*crossprod(Q*R,ylag)+phi)
    lambda <- betaS*iota[relevant]^alpha
    ll <- dnbinom(x=obs,mu=lambda,size=exp(-psi),log=TRUE)
    rv <- array(0,dim=dim(iota),dimnames=dimnames(iota))
    rv[relevant] <- ll
    rowSums(rv)/rowSums(rv!=0)
  }
}

likStouffer1 <- function (mle) {
  readRDS("rankmat1.rds") -> rr
  mle %$% {
    theta <- exp(theta)
    phi <- exp(phi)
    iota <- (N^tau1)*(theta*(rr^tau2)%*%ylag+phi)
    lambda <- betaS*iota[relevant]^alpha
    ll <- dnbinom(x=obs,mu=lambda,size=exp(-psi),log=TRUE)
    rv <- array(0,dim=dim(iota),dimnames=dimnames(iota))
    rv[relevant] <- ll
    rowSums(rv)/rowSums(rv!=0)
  }
}

likXRad1 <- function (mle) {
  readRDS("radmat1.rds") -> rr
  mle %$% {
    theta <- exp(theta)
    phi <- exp(phi)
    iota <- theta*(rr%*%ylag)+phi*N
    lambda <- betaS*iota[relevant]^alpha
    ll <- dnbinom(x=obs,mu=lambda,size=exp(-psi),log=TRUE)
    rv <- array(0,dim=dim(iota),dimnames=dimnames(iota))
    rv[relevant] <- ll
    rowSums(rv)/rowSums(rv!=0)
  }
}

poplim <- 1e5

data.frame(
  SV=likStouffer1(mle.stouffer1),
  CD=likCompDest(mle.compdest),
  G=likGravity(mle.grav),
  XR=likXRad1(mle.xrad1),
  N=N
) %>%
  name_rows() %>%
  dplyr::rename(town=.rownames) %>%
  dplyr::left_join(coords,by="town") %>%
  dplyr::filter(N<poplim) %>%
  tidyr::gather(model,value,-N,-town,-long,-lat) -> spatialLiks

expand.grid(m2=1:4,m1=1:4) %>%
  dplyr::filter(m1<m2) %>%
  dplyr::mutate(
           model1=c("SV","CD","G","XR")[m1],
           model2=c("SV","CD","G","XR")[m2]
         ) -> comps

The above plots are limited to cities of less than \(10^{5}\) inhabitants, since larger cities have few or no fadeouts.

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8    LC_NUMERIC=C            LC_TIME=en_GB.UTF-8    
##  [4] LC_COLLATE=C            LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
##  [7] LC_PAPER=en_US.UTF-8    LC_NAME=C               LC_ADDRESS=C           
## [10] LC_TELEPHONE=C          LC_MEASUREMENT=C        LC_IDENTIFICATION=C    
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] doMC_1.3.6       ncf_1.2-9        sp_1.4-2         scales_1.1.1    
##  [5] ggplot2_3.3.1    pomp_2.8.0.3     nloptr_1.2.2.1   bbmle_1.0.23.1  
##  [9] iterators_1.0.12 foreach_1.5.0    magrittr_1.5     reshape2_1.4.4  
## [13] plyr_1.8.6       knitr_1.28      
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4.6        highr_0.8           pillar_1.4.4       
##  [4] compiler_4.0.0      tools_4.0.0         digest_0.6.25      
##  [7] tibble_3.0.1        evaluate_0.14       lifecycle_0.2.0    
## [10] gtable_0.3.0        lattice_0.20-41     pkgconfig_2.0.3    
## [13] rlang_0.4.6         Matrix_1.2-18       yaml_2.2.1         
## [16] mvtnorm_1.1-0       xfun_0.14           coda_0.19-3        
## [19] withr_2.2.0         dplyr_1.0.0         stringr_1.4.0      
## [22] generics_0.0.2      vctrs_0.3.0         tidyselect_1.1.0   
## [25] deSolve_1.28        glue_1.4.1          R6_2.4.1           
## [28] bdsmatrix_1.3-4     rmarkdown_2.2       purrr_0.3.4        
## [31] ellipsis_0.3.1      codetools_0.2-16    htmltools_0.4.0    
## [34] MASS_7.3-51.6       colorspace_1.4-1    numDeriv_2016.8-1.1
## [37] stringi_1.4.6       munsell_0.5.0       crayon_1.3.4

References

Xia, Y., O. N. Bjørnstad, and B. T. Grenfell. 2004. Measles metapopulation dynamics: A gravity model for epidemiological coupling and dynamics. The American Naturalist 164:267–281.